function [ analysis_output,window_centers,window_counts,varargout ] = ...
    apply_length_resolved( ...
    input_bundle,keywords,L_min,L_max,windows,window_width,...
    complex_in,non_dropped, ...
    filament_function,window_function,varargin)
%APPLY_LENGTH_RESOLVED Applies analysis function in length-resolved manner
%   
%   ----------------------------------
%   % Regular length-resolved analysis (no bootstrapping)
%   [ analysis_output,window_centers,window_counts ] = ...
%   apply_length_resolved( ...
%   input_bundle,keywords,L_min,L_max,windows,window_width,...
%   complex_in,non_dropped, ...
%   filament_function,window_function)
%
%   Inputs
%
%   input_bundle - bundle from which result sections are read. The result
%   bundles can be passed either as a string (for a single result bundle),
%   or as a cell array with several strings (for more than one result bundle)
%
%   keywords - struct_array with fields any_keyword, with_keywords,
%   not_keywords. If empty [], no keywords are applied
%
%   complex_in - if true, complex results in length and velocity are
%   considered, otherwise they are not considered
%
%   non_dropped - is assigned true, all filaments are considered
%   irrespective of if they were removed using the Tracedropper or not. If
%   anything other than true, i.e. unassigned, empty [], or false, the
%   Tracedropper rejection is considered
%
%   L_min - start point for box car run
%   L_max - end point for box car run
%   windows - number of box car windows
%   window_width - width of box car window in fraction of the distance from
%   L_min to L_max
%
%   filament_function - function to get the necessary results from each
%   filament. If a cell array of functions is passed, several functions are
%   applied to Result Bundle
%
%   window_function - function to be applied to all results inside each
%   window
%
%   Optional Inputs
%
%   [ analysis_output,window_centers,window_counts ] = ...
%   apply_length_resolved( ...
%   input_bundle,keywords,L_min,L_max,windows,window_width,...
%   complex_in,non_dropped, ...
%   filament_function,window_function,'verbose')
%
%   Lets the routine print status reports to the command line
%
%   Outputs
%
%   analysis_output - cell array containing the analysis output. For each
%   length point of the running average procedure, a cell element is
%   allocated. This means that the resulting cell array will be of
%   dimension 1-by-windows, where windows is the number of windows in the
%   running average, as specified by the input argument windows
%
%   window_centers - numeric (scalar) array of dimension 1-by-windows, that
%   contains the center points of the length windows [micrometers]
%
%   window_counts - numeric (scalar integer) array of dimension
%   1-by-windows that contains the number of filaments in the length
%   windows
%
%   Example:
%
%   % Creation of a frame-to-frame velocity 2-dimensional histogram,
%   % resolved by length and velocity
% 
%   % Definition of keywords
%   % This assumes that a result bundle has been created from analysis
%   % results that are ordered into folders named 'day_1', 'day_2'
%   % etc. and inside these folders into folders for dilutions, one of them
%   % being names 'dil_1_4'.
%   keywords = struct;
%   keywords.any_keywords = {};
%   keywords.with_keywords = {'dil_1_4'};
%   keywords.not_keywords = {};
% 
%   % Define running average parameters
%   L_min = 0.4; %mum
%   L_max = 5.0; %mum
%   windows = 50; %number of running average windows
%   window_width = 0.05; % Fraction of the overall length range used as
%   %windows width
% 
%   complex_in = true; % Include complex results, important for short
%   %filaments
%   non_dropped = false; % Apply selection mask from TraceDropper procedure
% 
% 
%   % Definition of filament level function. This function will be applied to
%   % every section, and the results should arrive in either a cell array or a
%   % numeric array form. In this example, we get from each section a cell
%   % array with one element per filament in this section. Each element, in
%   % turn, contains all frame-to-frame velocities of an observed filament
%   filament_function = @(section) ...
%         {section.trace_results.frame_to_frame_velocities};
% 
%   % Definition of window level function. This function is used to determine
%   % the results for every running average window. The results will be stored
%   % as elements of a cell array, where the cell array has one element for
%   % each running average window. The results therefore are not restricted to
%   % scalar values, but can be any valid object.
%   % In this case we make a histogram of all frame-to-frame velocities from
%   % all filaments that are contained in the length window
% 
%   % Definitions to have the same histogram procedure in every length window
%   V_high = 7.0; hist_edges = linspace(0,V_high,16);
%   bin_centers = (hist_edges(1:end-1)+hist_edges(2:end))./2;
%   window_function = @(filament_results) ...
%         histc([filament_results{:}],hist_edges);
% 
%   % Choose an input bundle from disk
%   [load_file,load_path] = ...
%         uigetfile('*.mat','Select result bundle file from disk.');
%   result_bundle = [load_path load_file];
% 
%   [analysis_output,window_centers,window_counts] = ...
%         apply_length_resolved( ...
%         result_bundle,keywords,L_min,L_max,windows,window_width,...
%         complex_in,non_dropped, ...
%         filament_function,window_function,'verbose');
% 
%   % Plot the filament counts vs. running average window centers
%   subplot(1,2,1)
%   plot(window_centers,window_counts,'k.-')
%   xlabel('L[mum]')
%   ylabel('Filament count')
% 
%   % Compile 2D histogram image of frame-to-frame velocity vs. filament length
% 
%   image_array = zeros(numel(bin_centers),numel(window_centers));
%   for ll = 1:numel(window_centers)
%       % Access the histogram values for this running average window
%       current_bin_values = analysis_output{ll};
%     
%       % Remove the last bin, because histc happens to store junk in that bin
%       % by definition.
%       current_bin_values = current_bin_values(1:end-1);
%     
%       % Normalize with overall count
%       if sum(current_bin_values>0)
%           current_bin_values = ...
%               current_bin_values./sum(squeeze(current_bin_values));
%           image_array(:,ll) = current_bin_values;
%       end
%   end
% 
%   subplot(1,2,2)
%   imagesc(window_centers,[0 V_high],image_array)
%   set(gca,'YDir','normal')
%   colormap(gray)
%   xlabel('L[mum]')
%   ylabel('V_{f2f}')
%
%   ----------------
%   Optional output arguments for bootstrap confidence intervals
%
%   Bootstrap sampling only works for window_function functions that
%   return scalar values. If other results are returned, an error might
%   occur or erroneous results are returned!
%
%   [ analysis_output,window_centers,window_counts,bootstrap_CI ] = ...
%   apply_length_resolved( ...
%   input_bundle,keywords,L_min,L_max,windows,window_width,...
%   complex_in,non_dropped, ...
%   filament_function,window_function,varargin)
%
%   Example:
%
%   % Creation of mean velocity running average plot with 95% confidence
%   % interval
% 
%   % Definition of keywords
%   keywords = struct;
%   keywords.any_keywords = {};
%   keywords.with_keywords = {};
%   keywords.not_keywords = {};
% 
%   % Define running average parameters
%   L_min = 1.0; %mum
%   L_max = 2.0; %mum
%   windows = 50; %number of running average windows
%   window_width = 0.15; % Fraction of the overall length range used as
%   %windows width
% 
%   complex_in = true; % Include complex results, important for short
%   %filaments
%   non_dropped = false; % Apply selection mask from TraceDropper procedure
% 
%   % Definition of filament level function. This function will be applied to
%   % every section, and the results should arrive in either a cell array or a
%   % numeric array form. In this example, we get from each section a numeric
%   % array that contains the scalar trace-based velocity of each filament.
%   filament_function = @(section) [section.trace_results.trace_velocity];
% 
%   % Definition of window level function. This function is used to determine
%   % the results for every running average window. The results will be stored
%   % as elements of a cell array, where the cell array has one element for
%   % each running average window. The results therefore are not restricted to
%   % scalar values, but can be any valid object.
%   % In this case we make a histogram of all frame-to-frame velocities from
%   % all filaments that are contained in the length window
% 
%   % Definitions to have the same histogram procedure in every length window
%   window_function = @(filament_results) mean(filament_results);
% 
%   % Choose an input bundle from disk
%   [load_file,load_path] = ...
%     uigetfile('*.mat','Select result bundle file from disk.');
%   result_bundle = [load_path load_file];
% 
%   [analysis_output,window_centers,window_counts,CIntervals] = ...
%     apply_length_resolved( ...
%     result_bundle,keywords,L_min,L_max,windows,window_width,...
%     complex_in,non_dropped, ...
%     filament_function,window_function,'verbose');
% 
%   % Cast the analysis output cell array into a numeric array
%   velocities = [analysis_output{:}];
% 
%   % Plot the running average mean velocity
%   plot(window_centers,velocities,'k.-')
%   hold on
%   plot(window_centers,CIntervals(1,:))
%   plot(window_centers,CIntervals(2,:))
%   xlabel('L[mum]')
%   ylabel('Filament count')
%
%
%   ----------------
%   Optional output arguments for bootstrap confidence intervals and
%   bootstrap samples
%
%   Bootstrap sampling only works for window_function functions that
%   return scalar values. If other results are returned, an error might
%   occur or erroneous results are returned!
%
%   [ analysis_output,window_centers,window_counts,bootstrap_CI ] = ...
%   apply_length_resolved( ...
%   input_bundle,keywords,L_min,L_max,windows,window_width,...
%   complex_in,non_dropped, ...
%   filament_function,window_function,varargin)
%
%   Example:
%
%   % Creation of mean velocity running average plot with 95% confidence
%   % interval
% 
%   % Definition of keywords
%   keywords = struct;
%   keywords.any_keywords = {};
%   keywords.with_keywords = {};
%   keywords.not_keywords = {};
% 
%   % Define running average parameters
%   L_min = 1.0; %mum
%   L_max = 2.0; %mum
%   windows = 50; %number of running average windows
%   window_width = 0.15; % Fraction of the overall length range used as
%   %windows width
% 
%   complex_in = true; % Include complex results, important for short
%   %filaments
%   non_dropped = false; % Apply selection mask from TraceDropper procedure
% 
%   % Definition of filament level function. This function will be applied to
%   % every section, and the results should arrive in either a cell array or a
%   % numeric array form. In this example, we get from each section a numeric
%   % array that contains the scalar trace-based velocity of each filament.
%   filament_function = @(section) [section.trace_results.trace_velocity];
% 
%   % Definition of window level function. This function is used to determine
%   % the results for every running average window. The results will be stored
%   % as elements of a cell array, where the cell array has one element for
%   % each running average window. The results therefore are not restricted to
%   % scalar values, but can be any valid object.
%   % In this case we make a histogram of all frame-to-frame velocities from
%   % all filaments that are contained in the length window
% 
%   % Definitions to have the same histogram procedure in every length window
%   window_function = @(filament_results) mean(filament_results);
% 
%   % Choose an input bundle from disk
%   [load_file,load_path] = ...
%     uigetfile('*.mat','Select result bundle file from disk.');
%   result_bundle = [load_path load_file];
% 
%   [analysis_output,window_centers,window_counts,CIntervals,bs_samples] = ...
%     apply_length_resolved( ...
%     result_bundle,keywords,L_min,L_max,windows,window_width,...
%     complex_in,non_dropped, ...
%     filament_function,window_function,'verbose');
% 
%   % Cast the analysis output cell array into a numeric array
%   velocities = [analysis_output{:}];
% 
%   % Plot the running average mean velocity
%   plot(window_centers,velocities,'k.-')
%   hold on
% 
%   % Plot all bootstrap resamples of the length-resolved velocity mean
%   [number_of_samples,~] = size(bs_samples);
%   for bb = 1:number_of_samples
%       plot(window_centers,bs_samples(bb,:),'r--') 
%   end
%   xlabel('L[mum]')
%   ylabel('Filament count')



%% Check if verbose progress report is required

if ~isempty(varargin) && strcmp(varargin{1},'verbose')
    verbose = true;
else
    verbose = false;
end

%% Check if a selection function should be applied to eeach section

if nargin > 11 && ~isempty(varargin{2}) ...
        && strcmp(class(varargin{2}),'function_handle')
    % Apply a selection function to each section
    selection_function = varargin{2};
else
    selection_function = []; % Empty selection function leads to no selection
end


%% Parameters

nboot = 500; % Bootstrap samples

%% Preliminaries

% ---
% Input bundle acquisition

if verbose
    fprintf('Acquiring result bundles...\n')
end

if ~exist('input_bundle','var')
    input_bundle = [];
end

if isempty(input_bundle)
    % If an empty array is given as input_handle, choose an input bundle
    % from disk
    [load_file,load_path] = ...
        uigetfile(pwd,'Select result bundle file from disk.','*.mat');
    input_bundle = [load_path load_file];
    
elseif ischar(input_bundle)
    % If a string is given as the input bundle, check if the string points
    % towards a directory, and if so let user pick a result bundle file
    % starting from that directory. If the string points towards a
    % .mat-file, try to use this file as a result bundle file to load
    is_mat_file = exist(input_bundle,'file');
    
    if exist(input_bundle,'file') && is_mat_file==7 %means that path is a directory
        % Pick file dialog with pre-set start directory
        [load_file,load_path] = ...
            uigetfile(input_bundle,'Select result bundle file from disk.');
        input_bundle = [load_path load_file];
    end
        
end

% ---
% Make sure that keyword fields are functional

if ~exist('keywords','var')
    keywords = [];
end

if ~isstruct(keywords)
    keywords = struct;
    keywords.any_keywords = [];
    keywords.with_keywords = [];
    keywords.not_keywords = [];
else
    if ~isfield(keywords,'and_keywords')
        keywords.and_keywords = [];
    end
    if ~isfield(keywords,'with_keywords')
        keywords.with_keywords = [];
    end
    if ~isfield(keywords,'not_keywords')
        keywords.not_keywords = [];
    end
end

% ---
% Assign complex_in with false if it is anything but true
if ~exist('complex_in','var') || isempty(complex_in) || ~complex_in
    complex_in = false;
else
    complex_in = true;
end

% ---
% Assign parmaeters if they are not given as function inputs, and make sure
% they have functional values
if ~exist('L_min','var') || isempty(L_min) 
    L_min = 0;
end
if ~exist('L_max','var') || isempty(L_max)
    L_max = 14;
end
if L_min == L_max
    %if L_min is the same as L_max, the window has to be widened
    fprintf('L_min and L_max are the same, L_min lowered to 0.\n')
    L_min = 0;
end
if L_min > L_max
    % Inverse L_min and L_max so that L_min is smaller than L_max
    fprintf('L_max <L_min, reversed to fix it.\n')
    [L_max,L_min] = deal(L_min,L_max);
end
if ~exist('windows','var') || isempty(windows) || windows < 1
    windows = 8;
end
if ~exist('window_width','var') || isempty(window_width) || ...
        (window_width <= 0 || window_width >= 1)
    window_width = 0.2;
end


if ~exist('non_dropped','var') || isempty(non_dropped) || ...
        ~non_dropped
    non_dropped = false;
end

%% Check that functions were passed correctly

if ~iscell(filament_function)
    filament_function = {filament_function};
    filfun_iscell = false;
else
    filfun_iscell = true;
end
number_of_fil_functions = numel(filament_function);
if ~iscell(window_function)
    window_function = {window_function};
    winfun_iscell = false;
else
    winfun_iscell = true;
end
number_of_win_functions = numel(window_function);
if winfun_iscell ~= filfun_iscell
   error(['Filament and window level functions have been passed as' ...
       'different array types. Make both cell array or plain function handle.']) 
end
if ~winfun_iscell && ~filfun_iscell
    % If the filament as well as the window level function are both passed
    % as plain function handles, return the results not as a cell array,
    % but as a plain numeric array
    return_plain = true;
end
if number_of_fil_functions ~= number_of_win_functions
   error('Different number of filament and window level functions.') 
end



%% Get average filament lengths, trace velocities, filament level results

if verbose
    fprintf('Extracting filament velocities...\n')
end

% Pull trace velocities out of the queried result sections
pull_velocity = ...
    @(result_section) [result_section.trace_results.trace_velocity];
trace_velocity = extract_by_keywords(input_bundle, ...
    keywords.any_keywords,keywords.with_keywords,keywords.not_keywords, ...
    pull_velocity,non_dropped,selection_function);

if verbose
    fprintf('Extracting filament lengths...\n')
end

% Pull filament average lenghts out of the queried result sections
pull_average_length = ...
    @(result_section) ...
    [result_section.trace_results.average_filament_length];
average_filament_length = extract_by_keywords(input_bundle, ...
    keywords.any_keywords,keywords.with_keywords,keywords.not_keywords, ...
    pull_average_length,non_dropped,selection_function);

if verbose
    fprintf('Extracting filament results...\n')
end

% Apply the filament level analysis functions to get filament level results
filament_results = cell(1,number_of_fil_functions);
for ff = 1:number_of_fil_functions
    filament_results{ff} = extract_by_keywords(input_bundle, ...
        keywords.any_keywords,keywords.with_keywords,keywords.not_keywords, ...
        filament_function{ff},non_dropped,selection_function);
end

% Remove complex solutions if requested
number_of_sections = numel(average_filament_length);
for cc = 1:number_of_sections
    if ~complex_in
       inds = find(imag(average_filament_length{cc})==0 ...
           & imag(trace_velocity{cc})==0);
       average_filament_length{cc} = average_filament_length{cc}(inds);
       trace_velocity{cc} = trace_velocity{cc}(inds);
       for ff = 1:number_of_fil_functions
           filament_results{ff}{cc} = filament_results{ff}{cc}(inds);
       end
    end
end

% Pool results across all sections
pooled_trace_velocity = [trace_velocity{:}];
pooled_average_filament_length = [average_filament_length{:}];
pooled_filament_results = ...
    cellfun(@(results)[results{:}],filament_results,...
    'UniformOutput',false);

%% Execute the box averaging procedure

if verbose
    fprintf('Calculating length window level features...\n')
end

% ---
% Specify the window lower and upper bounds

interval_length = L_max-L_min; % Length of whole interval
absolute_window_width = interval_length.*window_width;
window_lower_bounds = ...
    linspace(L_min,L_max-absolute_window_width,windows);
window_upper_bounds = ...
    linspace(L_min+absolute_window_width,L_max,windows);
window_centers = ...
    window_lower_bounds + absolute_window_width./2;

% ---
% Execute window level function to calculate window results

% Allocate array to contain the analysis results
analysis_output = cell(1,number_of_win_functions);

% Allocate array to count the filaments in each window
window_counts = zeros(1,windows);
for ff = 1:number_of_win_functions
    analysis_output{ff} = cell(1,windows);
end

for ww = 1:windows
   window_inds = ...
       pooled_average_filament_length>=window_lower_bounds(ww) & ...
       pooled_average_filament_length<=window_upper_bounds(ww);
   for ff = 1:number_of_win_functions
       analysis_output{ff}{ww} = ...
           window_function{ff}(pooled_filament_results{ff}(window_inds));
   end
   window_counts(ww) = sum(window_inds);
end

if return_plain
    analysis_output = analysis_output{1};
end

%% Calculation of bootstrap confidence intervals for window function


if nargout > 3
   % If confidence intervals are requested as a function output, bootstrap
   % results by including different sections

   if verbose
       fprintf('Determining bootstrap confidence intervals...\n')
   end

   boot_function = @(indices) ...
       boot_boxcar_analysis(windows,window_lower_bounds, ...
       window_upper_bounds,window_function,...
       average_filament_length,filament_results,indices);
      
   if nargout > 4
       varargout = cell(1,2);
       [boot_ci,bootstats] = bootci(nboot,boot_function,...
           1:number_of_sections);
       varargout{1} = cell(1,number_of_win_functions);
       varargout{2} = cell(1,number_of_win_functions);
       for ff = 1:number_of_win_functions
           starting_ind = 1 + (ff-1).*windows;
           varargout{1}{ff} = ...
               boot_ci(:,starting_ind:starting_ind+windows-1);
           varargout{2}{ff} = ...
               bootstats(:,starting_ind:starting_ind+windows-1);
       end
       if return_plain
           varargout{1} = varargout{1}{1};
           varargout{2} = varargout{2}{1};
       end

   else
       [boot_ci] = bootci(nboot,boot_function,...
           1:number_of_sections);
       varargout{1} = cell(1,number_of_win_functions);
       for ff = 1:number_of_win_functions
           starting_ind = 1 + (ff-1).*windows;
           varargout{1}{ff} = ...
               boot_ci(:,starting_ind:starting_ind+windows-1);
       end
       if return_plain
           varargout{1} = varargout{1}{1};
       end

   end
   
end


function analysis_output = ...
    boot_boxcar_analysis(windows,window_lower_bounds, ...
    window_upper_bounds,window_function,...
    filament_length,filament_results,indices)


number_of_win_functions = numel(window_function);


% Pool data from sections according to indices passed by bootstrap
% procedure
filament_length = [filament_length{indices}];
for ff = 1:number_of_win_functions
    filament_results{ff} = [filament_results{ff}{indices}];
end


% Allocate array to contain the analysis results
analysis_output = zeros(windows,number_of_win_functions);
% Generate combined feature extraction function
window_function_implicit = @(filament_results,window_inds) ...
    combined_window_function(window_function,filament_results,window_inds);

for ww = 1:windows
    window_inds = ...
        filament_length>=window_lower_bounds(ww) & ...
        filament_length<=window_upper_bounds(ww);
    analysis_output(ww,:) = ...
        window_function_implicit(filament_results,window_inds);
end


function window_results = ...
    combined_window_function(window_function,filament_results,window_inds)

number_of_win_functions = numel(window_function);
window_results = zeros(1,number_of_win_functions);
for ff = 1:number_of_win_functions
    window_results(ff) = ...
        window_function{ff}(filament_results{ff}(window_inds));
end